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Iron sesquilayers grown at room temperature on W(110) 
exhibit a pronounced coercivity maximum near a coverage of 
1.5 atomic monolayers. On lattices which faithfully reproduce 
the morphology of the real films, a kinetic Ising model is uti- 
lized to simulate the domain-wall motion. Simulations reveal 
that the dynamics is dominated by the second-layer islands, 
which act as pinning centers. The simulated dependencies of 
the coercivity on the film coverage, as well as on the tempera- 
ture and the frequency of the applied field, are very similar to 
those measured in experiments. Unlike previous micromag- 
netic models, the presented approach provides insight into the 
dynamics of the domain-wall motion and clearly reveals the 
role of thermal fluctuations. 

PACS Number(s): 75.70.Ak, 75.40.Mg, 64.60. Qb, 05.50.+q. 



I. INTRODUCTION 

Recently, there has been much-interest in ultrathin 
iron films on W(110) substrates. EteI The present work is 
mainly concerned with the so-called sesquilayers, which 
are films with coverages between one and two atomic 
monolayers. When grown at room temperature, such 
structures consist of a nearly perfect monolayer with 
compact islands of the second layer on top HI The mechan- 
ical and magnetic ^properties are profoundly affected by 
this morphology.&u For example, the coercivity exhibits 
a pronounced maximum as a function of the sesquilayer 
coverage. In particular, around a coverage of 1.5 mono- 
layers (ML) the coercivity exceeds that of a monolayer 
or a doublelayer by more than an order of magnitude.cl 
Similarly, clastic properties show an unusual behavior 
in this thickness range The reportedly observed ab- 
sence of magnetic long-range order in a certain region 
of coverage was first-interpreted as a manifestation of a 
spin-glass-like phaseu However, the required presence of 
frustrated antiferromagnetic interactions was difficult to 



justify, and later an. explanation based on surface rough- 
ness was proposed £ru 

The aim of the present article is to study a model of 
this system by means of computer simulation. Our goal 
is twofold. First, we provide a semiquantitative expla- 
nation for the experimental coercivity observations. In 
spirit, this is similar to the micromagnetic models pro- 
posed in Refs. || and ||, but in contrast to the essentially 
zero-temperature micromagnetic arguments it gives an 
insight into the dynamics of the magnetization reversal 
at nonzero temperatures. Our model reproduces with- 
out additional assumptions ±he experimentally observed 
temperatureB and frequencyu dependence of the coerciv- 
ity. Second, we want to demonstrate the usefulness of 
the Monte Carlo approach to the simulation of real mag- 
netic systems. This is an application-oriented extension 
of our previous woijk-,on magnetization switching in ki- 
netic Ising models,EHl3 where one can find the theoretical 
background for the present article. Some preliminary re- 
sults of this study have been discussed in Ref. |ll|. 

The remainder of this paper is organized as follows. In 
the next section we describe how the dynamic magnetic 
properties of iron sesquilayers are reflected by our com- 
putational model. Section III is devoted to simulation 
results and their comparison with experimental data. A 
summary of our results, as well as a brief discussion of 
the implications of the proposed model for the frequency 
dependence of the coercivity are given in Sec. IV. 



II. COMPUTATIONAL MODEL OF AN 
ULTRATHIN IRON FILM 

From the theoretical point of view, an iron monolayer 
deposited on a W(110) substrate is a nearly ideal ferro- 
magnetic Ising-like system.El Although this only implies 
that the critical fluctuations in the vicinity of the Curie 
temperature can be described by the Ising model, one is 
tempted to speculate that it could be used also away from 
the phase transition. However, what is more important 
at lower temperatures is the magnetic anisotropy. The 



1 



magnetocrystalline anisotropy of the distorted bee iron 
structure on the tungsten substrate is almost a hundred 
times larger than in bulk iron.0 The strength oLthe re- 
sulting anisotropy field, on the order of a few T,el makes 
the Ising model a reasonable approximation for this sys- 
tem. The Ising Hamiltonian is given by 



H 



J s i s j ~ I^H 51 Si 



(1) 



Here, the two-state variables Si — ±1 at the nodes of a 
computational lattice represent local magnetic moments 
pointing in the positive or negative in-plane easy-axis 
[110] directions, respectively. The first sum runs over 
all nearest-neighbor pairs in the lattice. One and two 
layers of a square lattice are used to model mono- and 
double-layer regions of the film, respectively. Therefore, 
the number of nearest neighbors for a given site depends 
on its environment. The ferromagnetic spin-spin inter- 
action, J = 8.73 meV, is fixed such that the critical 
temperature, T c — 230 K, of a monolayers is correctly 
reproduced by the exact value for our square computa- 
tional lattice. The second term, with the summation over 
all sites, represents the interaction with the component 
H of the external magnetic field along the easy axis, and 
p, is the magnetic moment of the computational spin S{. 
The latter depends on the model lattice spacing a, which 
we choose to be a = 6 A. With this value of a, a com- 
putational spin Si represents a part of an iron monolayer 
containing approximately 5 atoms and carries a magnetic 
moment of [i = 11.43 hb- We use the bulk value of the 
iron atom magnetic moment. 

The morphology of the iron sesquilayers has been well 
studied by scanning tunneling microscopy (STM), and 
high-resolution images of their surfaces are available. 
Since the roughness of the surface plays a crucial role in 
the magnetic behavior, we have chosen to digitize STM 
pictures from Ref. [l] to generate our computational lat- 
tices. This has the advantage that the lattices faithfully 
reproduce the real films. Figure [l] shows the morphology 
of the iron sesquilayers. 

The dynamics we use in our simulations is the stan- 
dard Glauber Monte Carlo dynamics with local updates 
at randomly chosen sites, which, can be related to the 
quantum mechanical dynamics .113 To speed up the sim- 
ulations in weak external fields, we use a variant of the 
rejection- free Monte Carlo algorithm described in Ref. [H^. 
The simulation time is measured in Monte Carlo Steps 
per Spin (MCSS). Although its precise relation to phys- 
ical time is not known, it is expected that one MCSS 
roughly corresponds to a typical inverse phonon fre- 
quency. Here, we fix the simulation time scale such that 
1 MCSS = t MC = 10~ 12 s. 




FIG. 1. Illustration of the sesquilayer morphology at a 
coverage of 1.26 ML, with snapshots of the propagating do- 
main wall at the times marked by arrows in Fig. 2. Time 
increases from top to bottom in the picture. The patches are 
the islands of the second iron monolayer on top of the first 
one. The area shown is 1170 Ax 610 A, and the island con- 
figuration was digitized from Fig.l j) of Ref. jjj. The size of 
the computational lattice was 195 x 102 x 1(2) (two layers for 
islands). T = 1.3 J/k B « 132 K, = 0.02 J or H w 0.26 T. 
The high-contrast region is the stable magnetic phase. The 
double arrow shows the easy-axi s direction. An animation of 
this simulation can be found at tittp : //www. scri . f su. edu/j 
rikvold. 
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Before describing how we estimate the coercivity from 
simulations, let us consider the possible magnetization 
reversal modes in real systems. There are two extremes 
for magnetization switching: driven either by nucleation 
or by domain-wall motion. What is observed in experi- 
ments on these ultrathin iron films is the second mech- 
anism. On a microscopic scale, the sample is "infinite," 
and somewhere there exist seeds of the stable phase which 
start to grow almost instantaneously when the field is re- 
versed. As a result, magnetization reversal in a typical 
region of the sample is caused by a domain wall which 
propagates across the observed area. When the driving 
magnetic field is weak and/or the disorder is strong, the 
domain wall is pinned and does not move on the exper- 
imental time scale. Thus, the coercivity can be identi- 
fied with the field which makes the domain wall mobile 
under given circumstances (disorder, temperature, time 
scale, ...). To exclude the nucleation time from our sim- 
ulations, we prepare the initial state with all spins equal 
to +1, except for those in a narrow strip along one of the 
short sides of the lattice. These spins are initialized to 
— 1 and represent a domain of the stable magnetization 
that has just propagated into the region under observa- 
tion. We apply a negative magnetic field and measure 
the time needed, t sw , to reach the state with 90% of the 
spins reversed. Since the reversal occurs by propagation 
of the domain wall, the change in the the reversed volume 
divided by the lattice width gives us the average distance 
D the domain wall traveled. Together with the elapsed 
time it can be used to estimate the average domain-wall 
velocity. 



III. SIMULATION RESULTS 

Figure [j] illustrates the propagation of a domain wall 
through a sesquilayer system as it was recorded during 
a particular simulation run. One notices that the do- 
main wall typically does not cross the second-monolayer 
patches. Instead, it feels an island as an obstacle and 
gets pinned near its boundary. To overturn the spins in 
the island, it is necessary to overcome a free-energy bar- 
rier which depends on the driving field, the temperature, 
and also on the shape and size of the island. Due to 
thermal fluctuations, the domain wall tries to enter the 
islands, but is usually driven back by the free-energy bar- 
rier. A "successful" fluctuation must be large enough to 
have greater probability to grow through the island than 
to shrink and vanish; a fluctuation which just meets this 
criterion is called critical. 

To estimate the size of a critical fluctuation, denote by 
I the linear extent of the region which reverses its mag- 
netization during a jump of the domain wall to its next 
metastable configuration. The interface energy change 
SEs is proportional to £, while the "volume" energy 
change SEh is proportional to Hi 2 . While the latter 
is negative, the interface energy change can be negative 



or positive depending on whether the domain- wall length 
decreases or increases. The mobility of the domain wall 
is essentially given by the slower processes, for which 
5Er\ > and SEh < 0. The new configuration will be 
stable with respect to the previous one only if the volume 
energy change is sufficient to compensate for the interface 
energy change. From this it follows that I ~ 1/H. Since 
the interface and volume energies of the intermediate un- 
stable domain-wall configurations scale in the same way, 
the free-energy barrier AF to be overcome will be in- 
versely proportional to the driving field, AF = S/(/xff). 
The field-independent quantity S depends on the local 
environment, as well as on the temperature. In such a sit- 
uation, nucleation theory predicts that the—waiting time 
needed to observe a critical fluctuation isoEl 



H 



-K 



exp 



(2) 



where the exponent K is taken to be equal -to 3 as ex- 
pected for the two-dimensional Ising model, □ and fee is 
Boltzmann's constant. Once a critical fluctuation is cre- 
ated, the domain wall reorganizes quickly until it reaches 
a new metastable configuration of lower free energy than 
the previous one. In this way, propagation of a domain 
wall through the sesquilayer is explained as a series of 
thermally induced nucleation events followed by rapid 
magnetization changes. 
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FIG. 2. Magnetization per spin vs. time during the same 
simulation from which the snapshots in Fig. 1 were taken (at 
the times indicated by arrows). Most of the time, the magne- 
tization remains constant, apart from small fluctuations. The 
movement of the domain wall occurs in jumps. Each jump 
represents a nucleation event, after which the domain wall 
advances quickly to its new metastable configuration. 

Figure |^ shows how the intermittent domain- wall mo- 
tion is reflected in the magnetization per spin. The hor- 
izontal parts of the plot represent situations in which 
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the domain wall remains unchanged, apart from small 
fluctuations. The steps correspond to events when criti- 
cal fluctuations appear, thereby causing small regions of 
the sample to reverse their magnetization into the sta- 
ble phase. The average time between the magnetization 
jumps naturally depends on the size of the sample ob- 
served. For large regions it is expected to scale as the 
inverse of the length of the domain wall. Similarly, the 
size of an individual magnetization jump is proportional 
to the fraction of the observation area that is switched 
during that jump. Therefore, the magnetization curve 
should appear smooth for observation areas much larger 
than a typical island size. 

It is interesting to note that similar jumps in the mag- 
netization have beep-, observed experimentally in some 
thin magnetic films.E-3 Although in the present case the 
typical size of these Barkhausen jumps is about ten times 
smaller than in that study, it is possible that the effect 
could be observed using similar magneto-optical tech- 
niques. 

While traversing the sample, the domain wall encoun- 
ters barriers of different heights, and the total time 
needed to travel a certain distance is the sum of con- 
tributions similar to Eq. (Q), but with different values 
of S. As our simulations show, the residence time in 
different metastable domain-wall configurations can be 
very different, suggesting that there is a broad distribu- 
tion of free-energy barriers. Since we cannot determine 
that distribution, we assume it to be uniform between 
zero and an upper cutoff A. This and Eq. @ lead us 
to an approximate formula for the time elapsed during 
the domain- wall propagation through a given sample (we 
only keep the dominant terms), 
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Here we have added the first term, which corresponds to 
the time elapsed during the "free" domain-wall propa- 
gation, when the velocity is proportional to the driving 
magnetic field H . This term is usually negligible, except 
in systems with coverage close to two. The parameters 
A, B and A are fit to the simulation data as described 
below. 

The effective domain- wall velocity u e g = D/t sw (with 
D standing for the traveled distance, defined as the 
switched volume divided by the width of the lattice), 
obtained from the averaged switching times, is shown in 
Fig. H as a function of the driving field for several different 
coverages of the sesquilayer systems. A qualitatively sim- 
ilar domain- wall velocity dependence on-drijjing field was 
found in experiments on magnetic films ,E£r£3 in which the 
domain-wall motion is also believed to be controlled by 
thermal fluctuations. 

What is shown in Fig. || is an analogue of a depinning 
transition rounded by a finite temperature. It is to be 
stressed that in the present case there is no sharp transi- 
tion, even in an arbitrarily weak field. Strictly speaking, 



the average domain-wall velocity is always nonzero, but 
it quickly approaches zero on a scale of field strengths 
which, in agreement with the experimental findings, de- 
pend nonmonotonically on the coverage. 
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FIG. 3. Average domain-wall velocity as a function of 
the driving magnetic field for several different coverages at 
T — 184 K. As the coverage increases towards ~ 1.56 ML, 
the field near which the velocity practically vanishes increases 
(a). For coverages closer to 2 ML, it decreases again with 
increasing coverage (b). 

We measured the switching times for several 
sesquilayer samples with various coverages as functions 
of the magnetic field. Using Eq. (||), we estimated the 
fitting parameters for each coverage. The resulting val- 
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lies are summarized in Table I. We note that the behavior 
of the fitting parameters as functions of the coverage is 
quite reasonable. Also, the effective velocity of the "free 
domain-wall motion" Dfj,H/A is,similar to what has been 
measured in pure Ising systemsJlj 

We want to point out the difference between our 
Eqs. ([| [?]) and an approach commonly used in the ex- 
perimental literatureJl3'EEIin which one supposes that the 
free-energy barriers decrease linearly with H. That would 
be appropriate only if the size of the critical fluctuations 
were fixed. The scaling argument leading to Eq. (|^) im- 
plies that that is not so in the present case. However, it 
is extremely difficult to distinguish numerically between 
linear and inverse linear dependencies if one only has data 
from a relatively narrow range of field strengths. That 
is usually the case in experiments, and we face the same 
problem here. However, attempting to fit our simulation 
results using a linear field dependence, we obtain fitting 
parameters which behave rather erratically as functions 
of the coverage. We therefore believe that the 
model" used here is more appropriate in the range of 
H studied here, also from the numerical point of view. 

Another question concerning the functional form of the 
velocity-field relation used here, is whether it may be 
extrapolated to fields sufficiently weak that the velocity 
vanishes on the experimental time scale. We know that in 
a relatively strong field the size of the critical fluctuation 
I varies as I ~ 1/H. But is this still true when / ex- 
ceeds the length scales characteristic of the island struc- 
ture of the film? Or, equivalently, is the disorder strong 
enough to provide for a kind of "weak pinning" on large 
length scales, with groups of islands acting as "collec- 
tive" pinning centers? Both qualitative arguments and 
exploratory simulations of simplified configurations indi- 
cate that this is the case. However, the coercivities ex- 
tracted from our simulations approximately correspond 
to magnetic fields in which I approaches the island size. 
We cannot exclude that in even weaker fields the relation 
of the velocity to H is characterized by parameters differ- 
ent from those obtained in the relatively strong-field re- 
gion to which our simulations are confined. Without be- 
ing able to perform simulations in extremely weak fields, 
the only way to judge the appropriateness of our extrap- 
olation is the comparison of our final results with the ex- 
perimental data. The good, semiquantitative agreement 
is an indication that our expectation is, indeed, correct. 

Having an approximation formula for the domain-wall 
velocity v c g(H) as a function of the field, we can estimate 
the coercivity H c . Since coercivity is usually measured 
in a time-dependent field, we have to translate our static 
measurements into such a situation. We restrict ourselves 
to relatively slowly varying and weak fields in which the 
domain-wall motion is the main magnetization switch- 
ing mechanism. Suppose that the system is subjected 
to a magnetic field Hosm(ujt). The average domain- wall 
displacement at time t is then J v e g(H sm(u;t'))dt' . De- 
note by Ls the typical distance between seeds of the sta- 



ble phase. Then, the coercivity can be estimated as the 
value of the driving field reached when the domain walls 
have traveled that distance, 



arcsin(i/ c / Ho) 



v eS (H sin(<t>))d<t>Kl . (4) 



The velocity v e g is given in dimensionless units, and the 
conversion factor to ms _1 , a/tyic, is shown explicitly. 
Note that this equation is determined up to a factor of 
order unity. Since v e g depends strongly on H, this does 
not affect the estimated coercivity significantly. It is the 
order of magnitude of the combination a/(tMcLs^) and 
the value of the barrier A which determine the coercivity. 

There are two quantities in our model which we cannot 
determine precisely, namely the typical distance Ls be- 
tween the stable-domain seeds, and the factor tyic which 
relates the Monte Carlo and the physical time scales. 
They always enter our formulas as a product, so there 
is a single undetermined parameter in the theory which 
can be regarded as a fitting parameter fixing the coerciv- 
ity scale. We use the following values for tMC, Ls, and 
u: t MC = 10" 12 s, L s = 10" 5 m, and ui = 2tt s" 1 . While 
these are "reasonable" values, they were simply chosen 
such that the coercivity simulated at 0.8T C (184 K) ex- 
hibits a maximum value (as a function of coverage) which 
roughly corresponds to the experimental observation!] 

We used H = 0.5 T for the value of the driving-field 
amplitude, "Wihrffi is about the field strengths used in 
experiments.Btru When the coercivity H c is much smaller 
than Hq, it depends only weakly on the amplitude and 
frequency of the driving field in the combination loHq. 
On the other hand, one observes much stronger depen- 
dence on the amplitude when H c and Hq are comparable. 
In the extreme case when Hq is less than the coercivity, 
there is no solution to our Eq. (^). This situation cor- 
responds to the, experimentally observed collapse of the 
hysteresis loopJa 

In Fig. |i| we plot the coercivity H c as a function of 
the sesquilayer coverage. The coercivity first steeply in- 
creases to a maximum at a coverage of about 1.4. Then it 
decreases as the coverage approaches a full double layer. 
This is very similar to the experimentally observed be- 
havior of the coercivity. As the temperature decreases, 
the coercivity rapidly increases. Sander et al. reported 
a coercivity higher than 0.3 T at 140 Ka Our simula- 
tion in that temperature range nicely reproduces that 
observation. To illustrate the effect of the driving-field 
frequency uj, we also show a coercivity curve (filled cir- 
cles) predicted for a measurement at a frequency 1000 
times higher than the one shown by open circles. 

It is expected that another coercivity maximum would 
be observed for coverages between 2 and 3, although it 
should be much less pronounced. □ We do not have suit- 
able experimental data to create lattice structures for this 
coverage region, but simulations on structures which sim- 
ply repeat the morphology of the thinner films do suggest 
the existence of the second maximum. 
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FIG. 4. Coercivity, H c , of the iron sesquilayer as a func- 
tion of coverage. Different curves represent measurements at 
different temperatures and frequencies. The lower curve is to 
be compared with Fig. 3 of Ref. ^. Representative error bars 
are shown to give a rough idea of the statistical uncertainty 
due to errors in the fitting parameters. 



IV. CONCLUSION 

In conclusion, the good agreement between the shapes 
and the dependencies on T and uj of the simulated and ex- 
perimental plots of coercivity vs. coverage suggests that 
a simple kinetic Ising model can describe the basic pro- 
cesses involved in domain-wall motion at nonzero temper- 
atures in highly anisotropic films. Our simulations cor- 
roborate that the coercivity maximum observed for iron 
sesquilayers on W(110) is caused by domain- wall pinning 
at the islands of the second monolayer. Although the 
pinning mechanism here is in principle similar to the one 
proposed in the micromagnetic model of Ref. |], the cru- 
cial difference is in the role of thermal fluctuations, which 
the present model takes into account. Thermal fluctua- 
tions allow the domain wall to overcome obstacles cre- 
ated by the second-monolayer islands by creating critical 
fluctuations of the domain-wall shape. The domain-wall 
motion is intermittent with long periods of quiescence 
punctuated by nucleation events which are followed by 
rapid forward jumps of the wall. Without any further 
assumptions, our model also reproduces an increase in 
the coercivity at lower temperatures. It is to be stressed 
that there is only one unknown parameter in the model. 
Once it is fixed, simulations at different temperatures and 
film coverages can be used to predict the coercivity. 

Another interesting point is that the model also pre- 
dicts how the coercivity depends on the sweep rate (or 
frequency) of the magnetic field in which the measure- 
ment is done. By integrating Eq. (|j) in the low-frequency 
limit one can show that the inverse coercivity depends ap- 



proximately linearly on the logarithm of the sweep rate 
(frequency). Figure ||(a) shows that such an approxima- 
tion is quite reasonable. However, the same data plot- 
ted on a \og(H c ) vs. log(u;) scale in Fig. ||(b) also pro- 
vide a quite good linear plot. Recent experiments by 
Suen et al.H show that the logarithm of the hysteresis 
loop area (which, because of the shape of the loop, can 
be replaced by the coercivity in this case) increases lin- 
early with log(w) over several decades. Thus, our obser- 
vation is in good agreement with experimental findings, 
but it also demonstrates that it can be difficult to draw 
conclusions about the functional form of the frequency 
dependence, even with data spanning several frequency 
decades. 
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FIG. 5. The proposed model predicts that the coerciv- 
ity depends on the time scale of the measurement. Approxi- 
mately, the inverse of the coercivity is linear in the logarithm 
of the magnetic field sweep rate (frequency), as shown in panel 
(a). The same data is shown also as log(// c ) vs. log(w) in 
panel (b) . Note that the plots are nearly linear in both cases. 
Points represent calculations based on our simulation data, 
and straight lines are merely guides to the eye. With uio = 2tt, 
the shown frequency range roughly corresponds to the one in 
Ref. g. 
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Finally, we note that the proposed model is interest- 
ing not only from the point of view of its application to 
sesquilayer iron films. It describes the basic physics of 
thermally activated domain-wall motion, which is com- 
mon to many experimental systems (for example the per- 
pendicularly magnetized ultrathin films studied in Rcfs. [ 
[u} |l^,|l|,|^|) for which, when suitably modified, the ki- 
netic Ising model could provide a tool for a better under- 
standing of many effects. Ising-like models with a higher 
number of states can.be used to study rough films with 
a four-fold anisotropyEi or can serve as an alternative to 
the spin-1 mean-field-type description-of the influence of 
the island structure on the coercivity.E3 The dynamics of 
the model is also closely related to more general-topics 
such as interface propagation in disordered mediaE3 and 
depinning transitions in the presence of thermal fluctua- 
tions or hysteretic creep of elastic manifolds. Ell 
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coverage 


AID [J] 


B/D [J 2 ] 


A/k B T [J] 


0.95 




0.08(1) 


0.007(7) 


1.13 




0.08(1) 


0.033(2) 


1.26 




0.16(1) 


0.075(5) 


1.39 




0.15(3) 


0.100(15) 


1.56 




0.22(2) 


0.096(3) 


1.69 


5.1(5) 


0.10(1) 


0.090(8) 


1.84 


6.2(3) 


0.12(2) 


0.059(5) 


1.95 


4.5(4) 


0.12(2) 


0.015(4) 


2.13 


5.1(5) 


0.09(2) 


0.026(5) 



TABLE I. Fitting parameters of v e g = D/t sv , for several 
coverages at T = 184 K. Because of insufficient data acuracy, 
only a two-parameter fit was used for low coverages, where the 
term proportional to A is small compared to the exponential 
contribution in Eq. (^). Numerical uncertainities were esti- 
mated by processing two or more independent data sets. 



7 



